Numerical approach to low-doping regime of the t- J model 
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We develop an efficient numerical method for the description of a single-hole motion in the antiferromagnetic 
background. The method is free of finite-size effects and allows calculation of physical properties at an arbitrary 
wavevector. Methodical increase of the functional space leads to results that are valid in the thermodynamic 
limit. We found good agreement with cumulant expansion, exact- diagonalization approaches on finite lattices 
as well as self-consistent Born approximations. The method allows a straightforward addition of other inelastic 
degrees of freedom, such as lattice effects. Our results confirm the existence of a finite quasiparticle weight near 
the band minimum for a single hole and the existence of string-like peaks in the single-hole spectral function. 

PACS numbers: 71.10.Fd,71.10.Pm,74.25.Jb,79.60.-i 



INTRODUCTION 

A description of hole motion in the antiferromagnetic 
(AFM) background as described by the t — J model rep- 
resents one of the long-standing, open, theoretical problems 
in the field of correlated systems. The accurate solution of 
this problem may be crucial for understanding the behavior of 
high-temperature superconductors in the underdoped regime. 
Apart from the analytical solution in the Nagaoka regime UJ] 
in the limit of zero doping and small J/t, as well as rigor- 
ous theorems in the symmetric point J = 2t there are 
no exact solutions of this model in two spacial dimensions. 
Many outstanding, early approaches to this problem, such as 
the self-consistent Born approximation (SCBA) 
self-consistent perturbational approach (SCPA) j^, calcula- 
tions based on the string picture l3, 0, fioll . cumulant expan- 
sion (CE) technique 111 111 , exact diagonalization (ED) calcu- 
lations on small clusters \12, Tsi \]M , quantum Monte carlo 
calculatio ns [1 1511. and recent state-of-the-art QMC calcula- 
tions iflil IitL llSl] have provided quantitative description of 
the quasiparticle band- width, effective mass, and quasiparticle 
weight. Most of these methods reproduce dynamical proper- 
ties, such as the one-hole spectral function, as well. 

Among these approaches, ED calculations on small clusters 
provide exact solutions of the t-J Hamiltonian but suffer from 
finite-size effects due to small system sizes. Similarly, QMC 
calculations are limited to small, even though larger clusters. 
In addition, analytic continuation is necessary to obtain spec- 
tral properties, since most of QMC methods compute Green's 
function, defined in imaginary time. The SCBA and SCPA 
calculations are likewise limited to finite-size calculations in 
momentum space. Furthermore, they seem to overempha- 
size the string effect. On the other hand, early calculations 
based on the string picture are similar to the concept of the 
linear combination of atomic orbitals, which provide results 
for arbitrary momentum transfer. However, previous results 
HESH are not necessarily comparable with SCBA, SCPA 



and ED results of the t-J model. This is predominantly due 
to a limited number of variational parameters or to the small 
size of the Hilbert space used in these calculations. In com- 
parison to the variational approach used in Ref. lfl9l] . where 
the authors use a similar method for construction of the func- 
tional basis set, our method employs an exact-diagonalization 
approach using Lanczos technique, which allows solutions of 
much larger Hilbert spaces. 



The aim of this work is to present an accurate exact diag- 
onalization method, defined over a limited functional space 
(EDLFS). The method is based on the string picture Isjlgi llOll . 
which provides solutions to a single-hole problem in the AFM 
background that are free of finite-size effects. Furthermore, 
the method takes the advantage of modern computing capabil- 
ities that allow solutions of large matrices. Through the effi- 
cient construction of the limited functional space (LFS), even 
when using only a few thousand states, this method provides 
results that can be directly compared to the state-of-the-art nu- 
merical approaches on small lattices LI 41 that require tens of 
millions of states. 



Despite much work in this area !l20ll there remain many 
open questions concerning the physics of a doped AFM in 
the zero-doping limit. Current interest in this field is in 
part focused on the influence of the electron-phonon interac- 
tion on correlated motion of a hole in the AFM background 
[17, 21, 22, 23]. Another open question concerns the proper 
description of the difference between the hole- and electron- 
doped cuprates |24]. There is also a need for a method that 
would resolve the issue of a disappearing quasiparticle weight 
that was predicted in a thermodynamic limit of a doped AFM 
due to the phase string effect 12511 . 
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METHOD 

We start by writing the t-J model as 

H = Ht + H\\+ i?_L 
Ht ^ -tY^ c! ^Cj,, + h.c. 

^11 J^^i^! 

(iJ> 

ffx = J/2Y,StS7 + SrSf, 

(iJ> 



where q s = Ci_s(l — is a fermion operator, projected 

onto a space of no double occupancy. The sums run over 
the pairs of nearest neighbors as well as over the two spin 
orientations. Following in part works by Trugman 1911, In- 
oue and Maekawa [19] . and El Shawish and J.B. 12611 . we 
construct the LFS starting from a Neel state with one hole, 
(/)|^"'') = Ck|Neel) and proceed with generation of new states, 
by application of the kinetic part of the Hamiltonian Ht, i.e. 



(2) 



This procedure generates strings with maximum lengths given 
by N}i. While constructing the LFS, translation symmetry, 
generated by two minimal translations ri 2 = is 
taken into account. Due to exponentially rapid growth of the 
LFS with increasing Nh we introduce an additional parameter 
-^6 < Nh that restricts generation of long strings by imposing 
a condition under which all coordinates of spin-flips should 
satisfy ~ /^/| < ^b;/^ = {^iV} where h and / refer 
to electron and spin-flip indexes, respectively. Application of 
this condition improves the quality of the LFS by increasing 
the number of states containing spin-flips in the vicinity of the 
hole while keeping the total amount of states within compu- 
tationally accessible limits. The full Hamiltonian in Eq. |2]is 
then diagonalized within this LFS using the standard Lanczos 
procedure. We also note that our method, even though defined 
on the infinite lattice, is variational. Increasing the number of 
LFS systematically lowers energies of the zero- and single- 
hole states, i.e. and E^^\ 

While generation of single-hole states through application 
of only the kinetic part of the Hamiltonian seems a rather nat- 
ural choice for the construction of the single-hole wavefunc- 
tion, there remains a question of how to construct the LFS for 
the undoped case, i.e., the Heisenberg model. The solution of 
the latter seems necessary in order to compute spectral prop- 
erties of the one-hole system as well as its energy, relative to 
the undoped case. We next assume that the spacial extent of 
the disturbance of the spin background around the doped hole 
(in the literature also referred as a magnetic polaron) is finite. 
In this case it is not necessary to obtain the exact solution of 
the undoped system on the infinite 2D lattice for the correct 
description of the single hole properties. It is sufficient to find 
a solution of the Heisenberg model in the vicinity of the doped 



hole. We therefore construct the 0-hole LFS using the 1-hole 
LFS by simply filling the empty space with a spin. 
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Figure 1: (Color online) a) Single-hole energy -Bk=(ir/2.7r/2) vs. J/t. 
The full line represents EDLFS results of the present work that were 
obtained using Nst = 37402972 states. Dotted line represent CE re- 
sults from Ref. ITlll . The dot-dashed line, given by i5k={7r/2.ir/2) ~ 
J/t - 3.24 + 2.65(J/t)° '^^ is an interpolation based on ED re- 
sults obtained on 32 sites fT4ll . The double-dot-dashed line repre- 
sents interpolation £^=(71/2^2) = — 3.36 + 3.50(J/t)^/"^ based on 



The dashed line, 



k={7r/2,7r/2 



WMC method from Ref.^d 

\O.73 



-3.17 + 2.83(J/t)° " is a result of SCBA calculation, 
Squares represent QMC calculations from Ref. 1 18]. Note that -Bk/i 
results of ED, WMC and QMC calculations were shifted by — J/t 
due to a different deffinition of the t-J model, b) Quasiparticle 
weight ^k={7r/2,ir/2) VS. J/t. The full line represents EDLFS 
calculation using Nst = 37402972 states. The dot-dashed line 
represents ED results as summarized in the extrapolation given by 



^k={7r/2,ir/2) 



-0.136 + 0.664(J/t)" 



LMl, 



while the dashed 

line represents SCBA calculation .Zk=(,r/2,,r/2) = 0.63( J/t)"'®®^ 
from Ref. fS] (similar extrapolation was also obtained using SCR4 
in Ref. |6]). Diamonds represent WMC calculation from Ref. fl^ . 
and squares are QMC results from Ref. 1 18]. Note that WMC and 
QMC data were multiplied by a factor of 2 due to a different defini- 
tion of ^k used in Refs. QlQl]. c) The bandwidth W/t vs. J/t 
calculated with EDLFS (full line), ED (H (dot-dashed line), SCBA 
|5](dashed lines), QMC calculations from Ref. fTlll (squares), and 
W/t = 2t^/J (dotted line), proposed in Ref. f?]. d)Z^ vs. J/t 
(ten nearly overlapping curves) obtained by using 10 different sizes 
of LFS as explained in the text. 



STATIC PROPERTIES 

We now turn to the numerical results. The one-hole en- 
ergy, measured from the energy of the undoped system, _Ek ~ 
— is shown at the one-hole band-minimum k = 
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Figure 2: (Color online) a) Single-hole energy iJk and b) quasipar- 
ticle weight Z]^. Both quantities were calculated at fixed J/t = 0.3 
and 0.4 as noted in the figure. EDLFS results of the present work 
were obtained with Nst = 37402972 states (full lines), ED results 
(circles), calculated on a 32 site system, are from Ref. 1 14], QMC 
results (squares) are from Ref. IISI . while SCBA results (diamonds) 
in b) are from Ref. fS]. 



(7r/2,7r/2) in Fig. [T}i, along with CeJTiI]. ED Q, worm 
quantum Monte Carlo (WMC) IitII . Quantum Monte 
Carlo (QMC) calculations III], and SCBA Jl [H results. 
While our method is defined on the infinite system, the ab- 
solute values of E^^ and are ill-defined since they grow 
with the increasing number of basis states as the number 
of spin-flips, generated by the hole motion, increases. In 
contrast, i?k remains finite and well defined. Although our 
method can not be directly compared to the cumulant expan- 
sion technique 111 ill , we use some of the aspects of this tech- 
nique. Since we use in our method only the hopping part of 
the Hamiltonian, Eq. |2] to generate new states, all spin-flips 
are by construction limited to the vicinity of the hole. This 
by no means restricts the LPS only to connected strings. A 
propagating hole can also generate disconnected strings. In 
our approach, the precision of the description of the quantum 
spin fluctuations, caused by the presence of the hole, increases 
with decreasing distance from the hole. We can therefore ex- 
pect to achieve a thermodynamic limit as soon as the extent of 
the spin-flips in the LPS exceeds the size of the magnetic po- 
laron. In addition, we should stress that the zero-hole energy, 
i?"'*, per se has no physical meaning. It simply represents the 
solution of the Heisenberg model, defined on the zero-hole 
LPS, that is identical to the one-hole LPS with the exception 
of the additional spin located on the hole position. The high 
efficiency of our approach is reflected in good agreement of 
our results with CE method, jUIl], and SCBA approach IsIIItIi. 
Por comparison we also present results obtained with ED cal- 
culation on a 32 site system tl4.1 . WMC calculations [16], as 



well as with QMC calculations performed on much larger 
lattices (24x24) (see Pig.[I}i). In general, EDLPS, CE as well 
as SCBA methods give consistently lower values of the single- 
hole (polaron) energy in comparison to ED and QMC meth- 
ods. Here we stress, that the single-hole energy is extremely 
sensitive to the appropriate choice of the LPS for the 1 - as well 
as of the 0- hole space. Our results can be almost perfectly fit- 
ted with a form -E(7r/2.7r.2)/^ = = ao + bo{J/ty"> where 
parameters a, b and k = are listed in the first row of TableU 
We present the bandwidth W in Pig.[T]; along with SCBA 
Jstl and QMC results IllSi] , as well as with analytical predic- 
tion W = 2t^/J H, valid in the large J/t limit. We find 
good agreement with ED results in the physically most rele- 
vant regime J/t ^ 0.4. We note that in Ref. [ 14i] is defined 
?&W = E^Tj_^-f — £'(ir/2,7r/2)- In our approach due to broken 
translation symmetry the point k = (tt, tt) is folded onto the 
k = point. We thus believe that our definition of W is com- 
parable to the one in Ref. Q- QMC results from Ref. Ill] 
in contrast to EDLPS, ED and SCBA results predict slightly 
larger values of W/ 1. Note however larger error bars in QMC 
results around J/t < 0.6. 




(7i/2,7c/2; 



Figure 3: (Color online) a) Surface plot of at J/t = 0.3, 
computed on a mesh of 400 k— points. Variational space with 
Nst = 5213618 was used. Only 1/4 of the AFM BZ is shown, b) 
Contour plot of with equidistant contours while wavevectors are 
spanning the whole AFM BZ. The rest is the same as in a). 

Our calculation of £'k presented in Pig. |2^ reflects another 
important advantage of the present method over ED calcula- 
tions on limited system sizes. Note, however, that our calcu- 
lations are limited to the reduced APM Brillouin zone (BZ) 
because of broken translational symmetry. Defining the LPS 
on an infinite lattice allows calculation of physical properties 
at an arbitrary wavevector, limited to APM BZ. In Pig.|2^, we 
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present the dispersion relation i^k, calculated at J/t = 0.3 
and 0.4. Taking into account the fact that we are computing 
_Ek in absolute units (we used no additional energy shift), we 
find good agreement with ED results 1 14], in particular when 
comparing the k— dependence of the single-hole energy and 
the bandwidth. We find qualitative agreement also with QMC 
results from Ref. 1 18], calculated at J/t = 0.4. QMC method 
predicts larger bandwidth as also seen in Fig. [T];. To further 
quantify the efficiency of our method, we present our results 
for i?k at selected k— points in Table HIl These results were 
obtained using different numbers of states spanning the LFS 
at the physically relevant value J/t ~ 0.3. It is encouraging 
that reasonable results for the one-hole energy can be obtained 
from a LFS as small as Ngt = 1121. 

So far, we have shown that our method is successful in ob- 
taining the ground-state energy of the magnetic polaron, how- 
ever, the current scientific interest and existing open problems 
primarily concern dynamic properties of a doped hole. Before 
moving to dynamic response, we next present our results of a 
closely related quantity, i.e. the quasiparticle weight, Zk, vs. 
J/t, see Fig.[TJ). We define Zk by 



(*"''l4ck|«'0'') ' 



(3) 



where "^f]^^ 



(vj/C"!) represent the ground state of the system 
with either one or zero holes. Note that the ground state 
has k = 0. The agreement with the ED calculation 
is surprisingly good for J/t < 0.3. The WMC calculation 
from Ref. |16] and QMC calculations from Ref. yield 



slightly smaller values for Z^. It is noteworthy mentioning 
that the two different QMC methods yield consistent values 
of Zk (note nearly perfect overlap between the two meth- 
ods at J/t — 0.4). We have also tested the J/t ^ oo 
limit where exact result based on the spin-wave approxima- 
tion from Ref. |28] yields Z = 0.822. Our method gives 
Z = 0.926, which can be further compared with SCBA result 
calculation that gives Z = 1. 

We next briefly discuss possible sources of errors affecting 
results obtained by different approaches. In case of SCBA 
calculations, the error is due to the approximate nature of 
the calculation since only non-crossing diagrams are taken 
into account. ED calculations are limited to small lattice- 
sizes that may lead biased results due to finite-size effects. 
QMC simulation from Ref. |18], based on the loop-cluster 
Monte carlo method for the AFM state and the hole prop- 
agation within the fixed spin background, yields increasing 
larger error bars as one approaches the physically relevant 
regime J/t < 0.6, while WMC method |16] suffers from he 
minus-sign problem. EDLFS naturally depends on the choice 
of the LFS. Increasing the number of LFS should yield re- 
sults that are free of finite-size effects and valid in the ther- 
modynamic limit. Nevertheless, a systematic error may oc- 
cur due to a particular algorithm used to create different LFS, 
Eq. 121 To demonstrate the stability of our results against the 
choice of different LFS, as well as a rapid convergence of 





f^n hn 'Jn 




-3.37 2.86 0.62 




-3.39 4.50 0.76 


U)2 


-3.12 5.56 0.72 



Table I: Fitting parameters of the lowest peak positions in 

^(7r/2,7r/2)(t^)- Fits have the from tj„ = a„ + 6,i(J/t)^". 



our method for 0.02 ^ J/t < 1 with increasing Ngt, we 
present in Fig. ([TJl) nine nearly overlapping curves depicting 
Zk. The curves were calculated using different LFS's with: 
Nst ^ 7610,9786,43884,80108,218950,642406, 912478, 
3109626 and 5213618 as obtained using LFS generator, Eq.|2] 
with various values of Nh and Njy. The close agreement of 
values for Zk given in Table HIl represents additional qualita- 
tive demonstration of convergence in our calculation. Note 
that results are only weakly dependent upon the choice of pa- 
rameters Nil and Ni, that define the generating algorithm for 
LFS. 

In Fig I2J5 we present Zk along the special symmetry lines 
in the reduced AFM BZ. The agreement with the ED result is 
good for large values of k = |k| while the agreement with 
SCBA calculation is poorer. The discrepancy between our 
method and SCBA is similar over the whole AFM BZ since 
the SCBA does not suffer from finite-size effects. Most im- 
portantly, we find the value of Zk to be very small around the 
k = point (see also Table |ll]l, followed by a sharp increase 
with increasing k. These observations are consistent with the 
SCBA result. 

A surface plot of Zk in Fig.[3t that consist of 400 fc-points, 
calculated on a system with Ngt = 5213618 states, shows 
the power of our method. As expected from results, plotted 
in Fig. 12] Zk shows a pronounced minimum located at k = 
followed by a rapid increase with increasing distance from the 
k = point. In Fig.fSj) we show contour plot of Zk over the 
whole AFM BZ where contour lines, representing values of 
Zk, are uniformly spaced in the interval [0, 0.3]. 



SPECTRAL FUNCTIONS 

In Figs.|4^-c we plot the hole spectral function ^k(^^), cal- 
culated at J/t — 0.3 using three typical values of k. We define 

A]f^{uj) as 



El^-E'"^)), (4) 



where |5'k'*„) and i?k'n represent excited states and energies 
of the 1-hole system. In many respects, our results agree with 
the ED calculations of Leung and Gooding, Ref. |14]. The 
quasiparticle peak is well defined for wavevectors lying on the 
edge of the AFM BZ. In particular, for ki = {tt/2, 7r/2) the 
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Nt 


Nst 


(tt/2 tt/2) 


Ek 


fO 0) 


(tt/2 tt/2) 


Zk 

(tt 0) 


CO 0) 
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4 


1121 


-2.01925 


-1.95213 


-1.44065 


0.29253 


0.32780 


0.00002 
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7610 


-1.99475 


-1.92799 


-1.47960 


0.32617 


0.33895 


0.03093 
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8 


9786 


-1.99951 


-1.92888 


-1.47982 


0.32349 


0.33803 


0.03077 


10 


4 


43884 


-1.98751 


-1.92209 


-1.45354 


0.32529 


0.35097 


0.03359 


10 


8 


80108 


-2.00182 


-1.92305 


-1.45542 


0.32486 


0.34104 


0.03098 


12 


4 


218950 


-2.00272 


-1.93757 


-1.46192 


0.32790 


0.34895 


0.03535 


12 


8 


642406 


-2.01059 


-1.92709 


-1.43991 


0.32633 


0.34345 


0.03009 


14 


4 


912478 


-2.00024 


-1.93322 


-1.47915 


0.32902 


0.34942 


0.03907 


14 


8 


4992876 


-2.01830 


-1.93175 


-1.44255 


0.32805 


0.34314 


0.02809 


14 


12 


5225818 


-2.01831 


-1.93175 


-1.44255 


0.32804 


0.34314 


0.02809 


16 


8 


37402972 


-2.02175 


-1.93205 


-1.44112 


0.32939 


0.34324 


0.02713 



Table II: and Z^, calculated for J/t = 0.3 and different sizes of the LFS, generated using different values of Nh and A*';,. 



also be distinguished at k2. At yet higher frequencies there is 
another well defined peak. It is located at lu2 ^ —0.81 at ki 
with the following scaling W2 = 02 + 62 («//^)^^, Table|T] This 
fit is valid in the regime 0.1 < J/t < 1.0. ED results from 
Ref. flT] also display a well defined but less sharp structure at 
these frequencies. Moving towards k2, this peak looses some 
weight, however, it remains well defined. Spectrum at larger 
Lu is broad and mostly featureless. We note different scaling 
with J/t between the quasiparticle peak, located at luq and 
string-like peaks, positioned at uji and 0J2, see Table U 

At ks = 0, i^) displays a much smaller quasiparti- 
cle peak at w ^ E^^ than found in ED calculations. The 
broad, mostly incoherent part moves to lower frequencies 
and slightly shrinks. Nevertheless, the incoherent structure 
is broader than in ED calculations. 

In Figs |4^,d-e we plot (oj) for 4 different sizes of the 
LFS, ranging from Nst = 37402972 down to 7610. A filling 
up of the incoherent part of the spectrum in the large lu interval 
0<a;/<<4is the predominant effect of increasing Nst- All 
special features in the range oj seem to be well captured 
within the smallest size LFS. In addition, spectral functions 
for the largest two LFS's (see Figs|4^ and d) nearly overlap in 
the whole ut/t regime. 

In Fig|5] we present the evolution of A]f^{uj) for k moving 
from (7r/2,7r/2) towards (0,0). All curves were computed 
using Nst = 5213618. Dark-shaded areas are graphic rep- 
resentations of Zk. The evolution of spectral functions with 
increasing values of J/t, calculated at k = (7i'/2, 7r/2), is pre- 
sented in Fig.|6] The two lowest string-like peaks are denoted 
with Roman numerals. Scaling as discussed in the beginning 
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Figure 4: (Color online) a-c) Spectral functions Ak{io) for three typ- 
ical values of k, d-f) j4k(<jj) at k = {tt/2, tt/2) calculated using 
different sizes of LFS's as indicated in the insets. In all cases, we 
have used J/t = 0.3 and an artificial damping of e = 0.05. 



peak is located at = cjo = i^ki , see also TableHIl The quasi- 
particle peak is well defined also for k2 = (tt, 0). In contrast 
to ED results, we see a tiny peak at ki, located at ui ^ —1.58 
that scales with J/t as 0^1=01 + bi{J/t)'^^, Table |T] This 
fit is valid in the regime 0.3 < J/t < 1.0. This peak can 
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(7t/2, 



7C/2) 




Figure 5: Spectral functions Au_{u!) for wavevectors k = 
(7r/2,7r/2) through (0,0). In all cases we have used J/t = 0.3, 
Nst = 5213618 and artificial damping e = 0.1. Dark-shaded areas 
are proportional to the quasiparticle weight Z^. 



of this section of the quasiparticle peak, and two lowest string- 
like peak positions, ujq,oji, and uj2, with J/t can be qualita- 
tively followed. With increasing J/t, the broad continuum at 
high- w/t transforms into well defined peaks. 



CONCLUSIONS 

In conclusion, we have developed an efficient numerical ap- 
proach for calculating physical properties of a doped AFM 
insulator in the zero-doping limit. The presented method is 
highly efficient, free of finite-size effects, and it allows for 
computation of physical properties at an arbitrary wavevector. 
EDLFS obviously has a few shortcomings: a) the method is 
limited to calculations in the zero-doping limit, b) due to the 
broken symmetry of the starting wavefunction, calculations 
are Umited to the reduced AFM BZ, and c) results depend 
on the number of states Nst spanning the LFS. However, for 
most static as well as dynamic quantities convergence to the 
thermodynamic limit with increasing Ngt can be achieved. 

Using EDLFS, we have computed the quasiparticle energy, 
quasiparticle weight, and spectral functions and compared 
values to known and established analytical as well as numeri- 
cal results. We found the best agreement with CE and SCBA 
calculations for the single-hole energy while for the quasipar- 



Figure 6: Spectral functions ^k=(ir/2,7r/2) i^) for various values of 
J/t are shown. We have used Nst = 5213618 and artificial damping 
e = 0.1. Dark- shaded areas are proportional to the quasiparticle 
weight Zk=(^/2,^/2). 



tide weight at J/t < 0.3 best agreement was found with ED 
calculations obtained from the largest system of 32 sites. Our 
method with an already small number of LFS produces results 
for static quantities, such as the energy dispersion, the band- 
width and the quasiparticle weight, that are directly compa- 
rable to the state-of-the-art ED calculations on small lattices. 
Our simulations show that the quasiparticle weight around the 
band-minimum k = {tt/2,w/2) remains finite in the ther- 
modynamic limit. The quasiparticle peak is separated by a 
pseudo-gap from well defined string-like peaks. Comparing 
our results to ED calculations, we find a much smaller quasi- 
particle weight at the k = point. 

Our method can be easily extended to compute other static 
as well as dynamic quantities, e.g., various correlation func- 
tions in the vicinity of the doped hole and optical conductivity. 
Furthermore, it allows for the inclusion of additional higher- 
order terms in the Hamiltonian, such as the next-nearest 
neighbor hopping term that allows comparison of hole vs. 
electron doped AFM systems. The method can be easily ex- 
tended to computation of bound two-hole properties by adding 
another hole to the LFS. EDLFS can also be adopted to com- 
puting single-hole properties of the i — J model on the trian- 
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gular lattice. In this case a different, 120° ordered Neel state, 
should be used for the starting wavefunction, Ref. |29] Finally, 
by adopting the method of Ref. 1 3^ 31], the present approach 
offers a natural extension to computation of the Holstein-t- J 
model. 

One of the authors (J.B.) acknowledges the warm hospi- 
tality during his visit at the Institute for Materials Research, 
Tohoku University, Sendai. J.B. also acknowledges stimulat- 
ing discussions with A. Ramsak, I. Sega, and R Prelovsek, M. 
Stout for providing editorial suggestions, J. Vidmar for en- 
couragement, and the financial support of the Slovenian Re- 
search Agency under grant Pl-0044. S.M. and T.T. acknowl- 
edge the financial support of the Next Generation Super Com- 
puting Project of Nanoscience Program, CREST, and Grant- 
in-Aid for Scientific Research from MEXT. 



[1] Y. Nagaoka, Phys. Rev. 147, 392 (1966). 

[2] S. Sorella, Phys. Rev. B 53, 15119 (1996). 

[3] A. Ramsak and R Prelovsek, Phys. Rev. B 42, 10415 (1990). 

[4] S. Schmitt-Rink, C. M. Varma, and A. E. Ruckenstein, Phys. 
Rev. Lett. 60, 2793 (1988). 

[5] G. Martinez and P Horsch, Phys. Rev. B 44, 317 (1991). 

[6] Z. Liu and E. Manousakis, Phys. Rev. B 45, 2425 (1992). 

[7] A. Ramsak and P Horsch, Phys. Rev. B 57, 4308 (1998). 

[8] W. R Brinkman and T. M. Rice, Phys. Rev. B 2, 1324 (1970). 

[9] S. A. Trugman, Phys. Rev. B 37, 1597 (1988). 
[10] B. L Shraiman and E. D. Siggia, Phys. Rev. Lett. 60, 740 (1988). 
[11] P Prelovsek, I. Sega, and J. Bonca, Phys. Rev. B 42, 10706 
(1990). 

[12] E. Dagotto, R. Joynt, A. Moreo, S. Bacci, and E. Gagliano, 

Phys. Rev. B 41, 9049 (1990). 
[13] D. Poilblanc, T. Ziman, H. J. Schulz, and E. Dagotto, Phys. Rev. 



B 47, 14267 (1993). 
[14] P W. Leung and R. J. Gooding, Phys. Rev. B 52, R15711 
(1995). 

[15] S. Sorella, Phys. Rev. B 46, 11670 (1992). 

[16] A. S. Mishchenko, N. V. Prokofev, and B. V. Svistunov, Phys. 

Rev. B 64, 033101 (2001). 
[17] A. S. Mishchenko and N. Nagaosa, Phys. Rev. Lett. 93, 036402 

(2004). 

[18] M. Brunner, F. F. Assaad, and A. Muramatsu, Phys. Rev. B 62, 
15480 (2000). 

[19] J. Inoue and S. Maekawa, J. Phys. Soc. Jpn. 59, 2110 (1990). 
[20] E. Dagotto, Rev. Mod. Phys. 66, 763 (1994). 
[21] A. Ramsak, P Horsch, and P Fulde, Phys. Rev. B 46, 14305 
(1992). 

[22] K. M. Shen, R Ronning, D. H. Lu, W. S. Lee, N. J. C. Ingle, 
W. Meevasana, F. Baumberger, A. Damascelli, N. P. Armitage, 
L. L. Miller, et al., Phys. Rev. Lett. 93, 267002 (2004). 

[23] O. Rosch, O. Gunnarsson, X. J. Zhou, T. Yoshida, T. Sasagawa, 
A. Fujimori, Z. Hussain, Z.-X. Shen, and S. Uchida, Phys. Rev. 
Lett. 95, 227002 (2005). 

[24] T. Tohyama, Phys. Rev. B 70, 174517 (2004). 

[25] D. N. Sheng, Y. C. Chen, and Z. Y. Weng, Phys. Rev. Lett. 77, 
5102 (1996). 

[26] S. E. Shawish and J. Bonca, Phys. Rev. B 74, 174420 (2006). 
[27] K. J. von Szczepanski, P. Horsch, W. Stephan, and M. Ziegler, 

Phys. Rev. B 41, 2017 (1990). 
[28] A. G. Malshukov and G. D. Mahan, Phys. Rev. Lett. 68, 2200 

(1992). 

[29] A. E. Trumper, C. J. Gazza, and L. O. Manuel, Physical Re- 
view B (Condensed Matter and Materials Physics) 69, 184407 
(2004). 

[30] J. Bonca, S. A. Trugman, and I. Batistic, Phys. Rev. B 60, 1633 
(1999). 

[31] J. Bonca, S. Maekawa, T. Tohyama, and P. Prelovsek, Work in 
progress (2007). 



